clear
%% Problem
%   'EqualMinima','Himmelblau','Sixhump','ModifiedRastrigin','Vincent',
%   'Shubert','Composition','IncreasingMinima','Rastrigin','Schaffer'.
ProblemSet.FuncType = 'ModifiedRastrigin';
ProblemSet.dimension = 3;
ProblemSet.k = [3,3,3];
% ProblemSet.BaseFunc = 3;

[d, x_bound, ObjFunc, OptSet, precision_init, precision, ~, MaximalBudget, FileName]...
    = problem_setting(ProblemSet);

% MaximalBudget = 10000;
n_max_candidate = [5,6,7,8,9,10,11,12,15,20,25];
iepsilon = 4;
epsilon = [0.1,0.01,0.001,0.0001];
OptSet.epsilon = epsilon(iepsilon);

%% problem
Problem.Dimension = d;
Problem.Domain = reshape(x_bound,[1,d*2]);
Problem.Sampling = @(n, region)Sampling_BoxConstraint(n, region, ObjFunc);
NewRegionNum = 2;
MaxPartitionDepth_d = ceil(log((x_bound(2,:)-x_bound(1,:))./precision_init)./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d);
Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
%% Algorithm
AlgorithmM.QuantileLevel = 0.3;
AlgorithmM.NewBudget = 3;
AlgorithmM.StopCriteria = [1,MaximalBudget];
ClearRadius0 = (x_bound(2,:)-x_bound(1,:))./(NewRegionNum.^(MaxPartitionDepth_d));
ClearRadius = 2*min(ClearRadius0); % the maximum distance within a non-partitionable region
AlgorithmM.radius = ClearRadius;
AlgorithmM.local_search = @(starts,step,MaximalBudget)LocalSearch_test(ObjFunc,x_bound,starts,step,precision(iepsilon)/2,OptSet,MaximalBudget);
%% OPTIMIZATION
rep = 100;
OptObtained = zeros(rep,length(n_max_candidate));
TotalBudget1 = zeros(rep,length(n_max_candidate));
TotalBudget2 = zeros(rep,length(n_max_candidate));
OptNoFound = zeros(rep,length(n_max_candidate));
Times = zeros(rep,length(n_max_candidate));
OptRemain = cell(rep,length(n_max_candidate));
%%
for j = 1:length(n_max_candidate)
    AlgorithmM.MaxSampleSize = n_max_candidate(j);
    AlgorithmM.n0 = ceil(AlgorithmM.MaxSampleSize/3);
    for i=1:rep
        disp(['n_max: ',num2str(n_max_candidate(j)),', replication: ',num2str(i),'.'])
        tic
        [optima, ~, SampleSet, ~, ls_budget, OptSetRemain]...
            = prsmmo_ls_test( Problem, AlgorithmM, OptSet );
        Times(i,j) = toc;
        disp(['time used: ',num2str(Times(i,j)),'.'])
        OptObtained(i,j) = size(optima,1);
        TotalBudget2(i,j) = ls_budget;
        TotalBudget1(i,j) = size(SampleSet,1) - ls_budget;
        OptRemain{i,j} = OptSetRemain;
        OptNoFound(i,j) = size(OptSetRemain,1);
    end
end
%%
clear i j iepsilon OptSetRemain SampleSet ls_budget optima
save(FileName)
a=[mean(TotalBudget1)',mean(TotalBudget2)',mean(TotalBudget1+TotalBudget2)',var(TotalBudget1+TotalBudget2)'];
%% figures
% plot_id = find(sum(OptNoFound==0,1) >= 98);
plot_id = 1:11;
plot_value = TotalBudget1+TotalBudget2;
% plot_value = size(OptSet.sol,1)-OptNoFound;
figure()
hold on
plot_value = plot_value(:,plot_id);
mean_plot_value = mean(plot_value,1);

plot(n_max_candidate(plot_id),mean_plot_value./max(mean_plot_value),'-o')
plot(n_max_candidate(plot_id),(mean_plot_value+norminv(0.975,0,1)*sqrt(var(plot_value)/rep))./max(mean_plot_value),'--')
plot(n_max_candidate(plot_id),(mean_plot_value-norminv(0.975,0,1)*sqrt(var(plot_value)/rep))./max(mean_plot_value),'--')
hold off
legend([FileName,': ',num2str(max(mean_plot_value))])
xlabel('','Interpreter','latex','String','$n_{\max}$ value')
ylabel('The average budget size')
set(gca,'Fontname','Times New Roman','FontSize',16);
